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Abstract. The metrization of the space of neural responses is an ongoing re- 
search program seeking to find natural ways to describe, in geometrical terms, 
the sets of possible activities in the brain. One component of this program are 
the spike metrics, notions of distance between two spike trains recorded from 
a neuron. Alignment spike metrics work by identifying "equivalent" spikes in 
one train and the other. We present an alignment spike metric having C p un- 
derlying geometrical structure; the £2 version is Euclidean and is suitable for 
further embedding in Euclidean spaces by Multidimensional Scaling methods 
or related procedures. We show how to implement a fast algorithm for the 
computation of this metric based on bipartite graph matching theory. 



1. Introduction 

The analysis of neural signals seeks to "translate" any set of neural impulses into 
a language we understand [Rick e~et al. 1997] . But as the neural signals elicited by 
the same stimulus are never exactly alike, we need a quantitative means of deter- 
mining when two neural signals serve the same purpose. In what follows, we shall 
restrict our attention to signals from isolated cell electrophysiology, and represent a 
"spike train" neural signal as a marked point process, an ordered list of spike times 
labeled by their respective neurons of origin |Aronov, Victor, 2004| . To compare 
spike trains, Victor and Purpura introduced the notion of a spike metric, a distance 
function on the set of spike trains endowing it with the topological properties of a 
metric space |Victor, Purpura, 1996| , |Victor, Purpura, 1997| . Spike metrics have 
been successfully used to quantify variability in data and characterize neural coding 
in the visual, auditory, olfactory, taste, and electric sensory systems |Victor, 2005| . 
Numerous spike metrics have been proposed with the goals of clustering neural 
signals by stimulus |Victor, Purpura, 1997 , [Schrauwcn, Campen hout, 2007] and 



embedding them in Euclidean space Victor, Purpura, 1997| via multidimensional 



scaling (MDS) Kruskal ' Wish, 1978| . MDS allows one to visualize the geometry of 
spike data, especially after the use of sophisticated dimensionality-reduction proce- 
dures such as Local Linear Embedding Roweis, S aul, 2000 . It also permits us to 



use more sophisticated clustering techniques designed to work only in vector spaces, 
such as PCA-guided K-means |Ding, He, 2004 and soft-margin support vector ma- 



chines, as done (without embedding) in Schrauwcn, Cam penhout, 2007] 



One approach to spike metric design for spike trains generated by single neurons 
uses the rate-coding hypothesis of spike generation, the idea that spike trains in- 
duced by the same stimulus are instances of the same variable-rate Poisson process 
[Rieke et al. 1997] , The implication is that to compare two spike trains, one must 
compare their estimated underlying rate functions. To estimate the rate function of 
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a spike train, represent it as a sum of Kroncckcr delta functions at the spike times 
and stream it through a low-pass filter, or, equivalently, convolve it with a nonneg- 
ative function, or kernel, that integrates to 1. Common choices of kernels include 
the one-sided exponential tail |van Rossum, 2 001 , the two sided exponential tail, 
the 'tent' function (0 above and below a certain range, and linearly increasing and 
then decreasing at the same rate within that range), and the normal distribution 
function (Schrauwen, Campenhout, 2007 . The value of the metric between two 
spike trains is the C p norm of the difference between their estimated rate functions 
(p> 1). For p = 2, such metrics resemble Euclidean distances, and MDS can embed 
them in (mostly) Euclidean finite dimensional vector spaces. 

However, their is substantial evidence that the time-coding of neural signals 
contains meaning beyond that conveyed by the firing rate |Chase, Young, 2006 1, 



jGerstner et al., 1997| , [Rieke et al. 1997] . |Samonds, Bonds, 2003 . Based on this 



fact, Victor and Purpura introduced a spike metric equal to the minimum 'cost' of 
aligning the spikes in two spike trains, |Victor, Purpura, 1996| , |Victor, Purpura, 1997| . 
Their approach was to generalize methods used for comparing strands of DNA. In 
particular, they reengineered Sellers' dynamic programming algorithm for align- 
ing and calculating the evolutionary distance between pairs of DNA sequences 
Sellers, 1974 to compute a distance between pairs of spike trains |Victor, Purpura, 1996| , 



Victor, Purpura, 1997| . This metric preserves the integrity of individual spikes 



instead of viewing them as contributions to a rate function. It has the added ad- 
vantage of being able to compute distances between spike trains that may have 
contributions from multiple neurons, and in the future, it may be generalized 
to align spikes from multiple spike trains at once, as has been done with DNA 
Notrcda me, 2002] , The problem with it is that it resembles an C\ norm on a vec- 
tor space. If one uses it to compute all the distances among a group of spike trains 
elicited from a common stimulus, and then embeds those spike trains using those 
distances by MDS, the resulting picture is very complicated. The embedded set 
of spike trains may have hyperbolic structure that is not present in the stimulus 
space. 

We propose a spike metric consistent with the time-coding hypothesis of spike 
generation that has all of the desirable properties of an C p norm. When p = 1, this 
metric is equal to the Victor-Purpura metric [Victor, Purpura, 1996| , |Victor, Purpura, 1997| . 
When p = 2, embeddings of spike trains using this metric 'fit' in Euclidean space 
with substantially less difficulty than they do for any other value of p > 1, so 
advanced dimensionality-reduction and clustering techniques may be used. The 
dynamic programming algorithm described by Victor and Purpura does not work 
when p > 1, so we propose a faster quadratic-time algorithm that works for all 
P > 1. 

This procedure we use to calculate our metric is the Hungarian Algorithm 
|Schrijver, 2003| , which can be made substantially faster for the specific task of 
comparing spike trains. The crucial insight behind the Hungarian Algorithm was 
independently discovered by Kuhn and Munkres |Kuhn, 1955 , [Munkrcs, 1957| . Its 



function is to solve the minimum-weight matching problem on weighted bipartite 
graphs. The Hungarian Algorithm is a special case of algorithms to solve gen- 
eral matching problems [Papadimitriou, Stciglitz, 1998] and assignment problems 
Burkard, 1999 1. Our version of the algorithm it is a special case of algorithms to 



solve transportation problems on Monge arrays, described in Burkard, Klinz, Rudolf, 1995| . 
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The original transportation problem was solved by Hoffman Hoffman, f963] and 
formulated by Monge [Mo nge, 1781| . 

2. The Metric 

We desire to match spikes in one spike train to corresponding spikes in another 
spike train to compare the two signals. Two trains fired in response to the same 
stimulus may have different numbers of spikes. So, an alignment of the spikes in 
the two trains will be inherently imperfect, some spikes may have to be deleted 
from the trains to make the alignment exact. Victor and Purpura's idea was to 
break the process of aligning two spike trains into steps, each with their own cost. 
The metric is equal to the sum of the costs incurred by the most efficient alignment 
|Victor, Purpura, 1996| , |Victor, Purpura, 1997| . Let q > 0. The cost of aligning 
two spikes on different trains is qAt, where At is the distance between them. The 
cost of deleting a spike on either train is 1. If in two trains the differences between 
spike times in aligned pairs of spikes are Ati, for i G {1, . . . ,fc}, the number of 
spikes deleted from the first train is D\ , and the number of spikes deleted from the 
second train is D 2 , the distance between the two spike trains is 

k 

^2qAt l + D 1 + D 2 . 

i=l 

The metric is defined to be the minimum of this quantity over all ways of aligning 
the trains. Our generalization is as follows: Let p > 1 , the metric is 

i/p 

5^ q p At P + D x + D 2 
_i=i 

When p = 1 this is the standard Victor-Purpura metric. 

We can formulate the metric in a way that is more mathematically specific as 
follows: Let x = = {xi} n L 1 and y = {j/j}™ =1 be spike trains, or strictly increasing 
finite sequences of real numbers. Let M be a 'matching' between the two spike 
trains, a set of ordered pairs {(x^ , yj 1 ), . . . , (xi k ,yj k )} with no element of x or y 
repeated. Let M. be the set of all matchings. Let the cardinality of M, \M\, be the 
number of ordered pairs contained in M. Define 

i/p 

]T qV\xi- yj \v + {m-\M\) + {n-\M\) 



(1) d p , q [M](x,y) 
Our metric is 



d p , q {x,y) = min . d p , q [M](x,y). 



We shall show below that we can compute this expression through the Hungarian 
Algorithm |Schrijver, 2003 , with additional speedups, from o(mn(m + n) 2 ) time to 



o(mn) time. 

3. Outline of the algorithm 

Our metric is a minimization over all possible matchings. We can map this 
minimization to a classical problem, weighted bipartite matching . In this problem 
we have a graph whose nodes belong to two different classes with no edge connecting 
members of the same class (bipartite) and each edge of the graph has a "weight" 
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or cost. A matching or pairing is a subgraph in which each node has at most one 
connection, and the matching with the smallest total cost is sought. 

A priori, the number of distinct matchings is enormous. There is a profound 
result, called the Hungarian Algorithm |Sc hrijver, 2003] , that shows that only a 
small subset of these matchings needs to be examined. This result works by showing 
that the optimal match with precisely k + 1 pairings (cardinality k + 1) is related in 
a rather specific way to the optimal match of cardinality k: the optimal k+1 match 
differs from the optimal k match by an "M-augmcnting path" . Since there are many 
fewer such paths than matchings of cardinality k+1 the search is considerably faster. 

We shall show that the convexity of our metric translates into a specific property 
of the graph weights called the Monge property Burkard, Klinz, Rudolf, 1995 . Us- 
ing the Monge property we can discard large fractions of the M-augmcnting paths 
as candidates: the paths have to be "monotonic" and cannot "jump" . We call the 
resulting, much smaller set of "M-augmenting incompressible monotonic paths" the 
shifts. There are far fewer shifts to explore in connecting level k with level k+1, 
at most m + n — 2k — 1 . 



4. The Shift Algorithm 

For certain matchings M we can construct a unique two-row matrix represen- 
tation, or 'matching matrix' as follows: The first row contains the elements of x 
in increasing order (from left to right), and the second row contains the elements 
of y in increasing order. If {xi,yj) € M, then Xi and yj are in the same column 
of the matrix. If Xi is not matched in M, then it is in the same column as a V 
symbol. The same goes for unmatched elements of y. Any unmatched elements of 
x and y together appear in increasing order from left to right in M. Making this 
matrix takes o(m + n) time. For instance, let x = {b,e,f, h}, y = {a, c, d, <?, i}, 
where a<b<---<h<i, with M = 0. The matching matrix looks like: 

b * * e f * h 
¥ c d * * g * 

Now assume '6' and l K' from x are matched with 'c' and 'g' from y. The new matrix 
looks like: 

* b * e / h * 
a c d * * g i 
Define the 'shift' operation on a matching matrix as follows: Delete two *'s, one 
on the top and one on the bottom and such that there exist no other *'s on either 
row between them. Then shorten the rows to compensate. Below are both possible 
shifts applied to the matrix above: 



/ h 



f 
fJ 



Main Theorem. If M minimizes d Pyq [M](x, y) over matchings of cardinality k < 
min(m, n), then at least one matching that minimizes d p q [M](x, y) over matchings 
of cardinality k+1 is some shift of the matching matrix of M. 

The proof is outlined in Section 6 and completed in the Appendix. The empty 
matching always has a matching matrix, and any shift applied to a matching ma- 
trix returns a matching matrix. Therefore, we can find the optimal matching of 
cardinality k + 1 by tabulating all shifts on the optimal matching of cardinality 
k and then picking the one that minimizes d p . (? [Af](x, y). Even if there are two 
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optimal matchings, M\ and Mi of cardinality k + 1, and only Mi is a shift of the 
matrix corresponding to the optimal matching of cardinality k, at least one optimal 
matching of cardinality k + 2 is a shift of Mi . 

The time cost of the naive implementation of the algorithm is as follows: Mak- 
ing the matching matrix takes o(m + n) time. Finding the distance-minimizing 
matching of order k + 1 by calculating the expected gain or loss to the distance 
function from every possible shift on the distance-minimizing matching of order 
k takes o(m + n) time. There are o(min(m, n)) rounds of searching, so we get 
o((m + n) min(m, n)) time. 

We have implemented two key speed-ups. min(m, n) rounds of shifting is the 
worst-case scenario; by Theorem 3 the algorithm can terminate as soon as all pos- 
sible shifts increase the value of d p q [M](x.,y). After making the first matching 
matrix, cut it in two between any two spikes (on either train) separated by more 
than 2 1 / p /q. By Corollary 4, the optimal matching is the union of the two opti- 
mal matchings found by repeatedly shifting these two submatrices. To prove the 
algorithm's correctness we need a result from graph theory. 

5. The Hungarian Algorithm 

For the sake of completeness we condense Chapter 3 of [Schrijvcr , 2003| below. 
Here we demonstrate one way in which the optimal matching of cardinality k + 1 
is related to the optimal matching of cardinality k. This is the first step toward 
showing that the relationship is in fact a shift. 

Definitions. A graph G = (V, E) is a set of vertices V and edges E, where an edge 
e is defined to be a pair of two vertices, e = {v,v'}, where v, v' G V. G is complete 
bipartite if V = A U B where A n B = and every edge e contains exactly one 
element of each of A and B. G is edge- weighted if there exists a weight function 
w : E — ► K on the edges of G. From now on all graphs mentioned will be complete 
bipartite and edge-weighted. A matching M on G is a subset of E such that no two 
edges in M share a vertex. Let w(M) = J2 eeM w{e) be the weight of matching M. 

A path P on G is an ordered subset of distinct elements of E, (ei, . . . , e t ), such 
that ei and e^+i share a vertex for all 1 < i < t but ei and et do not share a vertex. 
A path P is M-augmenting for a matching M if: 1. t is odd. 2. e2, e^, . . . , &t-\ <= M. 
3. ei and et each only share one vertex with an element of M. For an example of 
an M-augmenting path, see Figure 1. For an edge e, define its length function /(e) 
to be w(e) if e G M and -w(e) if e £ M. For a path P, let l(P) = J2 e eP l ( e )- 

If C,D C E, CAD contains all edges that are in exactly one of C and D. It 
follows that if P is Af-augmenting then MAP is a matching of cardinality |A/| + 1. 

Theorem 1. Schrijv er, 2003 1 , Section 3.5, Proposition 1: Let P be an M-augmenting 
path of maximum length. If M is a matching of minimum weight of cardinality k, 
then M' = MAP is a matching of minimum weight of cardinality k + 1 . 

The proof is in the appendix. 

6. Proof of Shift Algorithm 

Here we outline the proof that augmenting paths of maximum length correspond 
to shifts. The details of the proof are in the Appendix. 



6 ALEXANDER J. DUBBS 12 , BRAD A. SEILER 12 , AND MARCELO O. MAGNASCO 1 

Lemma 1. Computing d p>q (x, y) is equivalent to finding a minimum weight match- 
ing on the complete bipartite weighted graph with vertex sets A = {ai, . . . , a m } and 
B = {b\, . . . , b n } and weight function w(ai,bj) = q p \xi — yj\ p — 2. 

Proof. Every matching M on this graph is assocated with a matching M (abusing 
notation) between spike trains x and y. To minimize d p . q [M](x,y) over M it is 
equivalent to minimize 

]T (q p \xi- yj \ p -2)= Yl w{a i ,b j )=w{M) 

(li,Sj)£M (ai,bj)eM 

□ 

Now we need to find an analogue for shifting a matching matrix in terms of 
performing an operation on graphs. It turns out that the analogue is taking the 
symmetric difference between M and a special type of M-augmenting path, which 
we will also call a "shift" . It is a path of the form: 

{{<H,bj}, {ai, bj+i}, {fli+i, b j+ i} {a i+N , b jj+N }) 

or 

({ai, bj}, {a i+ i,bj}, {a i+1 ,b j+ i}, . . . , {a i+N , b j+N }), 

with a couple additional constraints: in the first case there are no unmatched 
elements of x between Xi and yj and there are no unmatched elements of y between 
X i+N an d Ui+N, an d in the second case there arc no unmatched elements of y 
between Xi and yj and there are no unmatched elements of x between Xi + N and 
yi+N- See Figure 2 for an M-augmenting path that is a shift. In this figure and 
all subsequent figures we place the sets of nodes corresponding to the spikes on 
trains x and y on parallel lines, and nodes farther to the right correspond to spikes 
occurring later in time than the spikes associated to nodes on their left. We can 
now restate the Main Theorem: 

Main Theorem, Restated. Let a bipartite graph have vertex sets A = {ai, . . . , a m } 
and B = {bi, . . . , b n } and weight function w(ai, bj) — q p \xi — yj\ p — 2. Let M min- 
imize w(M) over matchings of cardinality k < min(m, n). Then there exists an 
M-augmenting path P that is a shift such that MAP minimizes w(M) over match- 
ings of cardinality k + 1 . 

Lemma 2. Strict Monge Property: Let x il < Xi 2 and yj 1 < yj 2 . For p > 1, 
\ x h - Vn \ P + 1^2 - Vn \ P < \ x h - Vh \ P + \xi 2 - V n \ P 

The proof uses convexity and is in the Appendix. 

Lemma 3. Let p > 1. If M is a minimum weight matching of cardinality k 
containing edges {a il ,bj 1 } and {a i2 ,bj 2 } then either ii < i 2 and ji < j 2 or i\ > i 2 
and ji > j 2 ■ 

Proof. Remove edges {a^ , bj 1 } and {ai 2 , bj 2 } and replace them with edges {a^ , bj 2 } 
and {ai 2 , bj t } to form the matching M' of cardinality k. By Lemma 2, 

w(M') - w(M) = qP\x tl - y j2 \P + qP\x t2 - y h \ p - q p \x h - y n \ p - q p \x t2 - y 32 \ p 

Since M is of minimum weight, this expression must be positive. Lemma 2 implies 
that either ii < i 2 and ji < j 2 or i\ > i 2 and ji > j 2 . □ 
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This lemma can be pictorially understood as meaning if M is a minimum weight 
matching of cardinality k, it contains no edges that cross, using our convention for 
drawing figures stated below Lemma 1 . See Figure 3 for a demonstration of Lemma 
3. 

Theorem 2. Let p > 1, let M be a minimum weight matching of cardinality k, 
and let P be an M-augmenting path of maximum length. Then P is of the form 

i h h L . h n }> i a *2 . bj 2 }, • • • , , b jt }) 

or 

({Oil i & il }) { fl »2 . }> {^2 . 6 J2 }> • • • . l a *, . fy ( }) 

where the sequences {ai t } and {bi t } are either strictly increasing or strictly decreas- 
ing. 

The exact proof is in the Appendix. It can be understood using diagrams. What 
we want to show is that P does not turn back on itself, it either moves forwards 
or backwards through the graph. If it should decide to change direction, one of 
two things happen. It might make use of two edges in M that cross, contradicting 
Lemma 3, sec Figure 4a. If not, it uses two edges that are both not in M that 
cross (if it uses an edge in M and an edge not in M that cross, one of the previous 
two problems necessarily occur). In the second case, both of those edges will be 
in MAP. But if P is of maximum length, by Theorem 1, MAP is a minimum 
matching of cardinality k + 1, which also cannot contain crossing edges by Lemma 
3, see Figure 4b. This is a contradiction, and it follows that P cannot turn around. 

We are now ready to sketch the proof of the Main Theorem for p > 1. It follows 
for p — 1 by real analysis. Both proofs are in the Appendix. 

Let M be a minimum weight maching of cardinality k, and P be an M- augmenting 
path of maximum length. P is of the form described by Theorem 2. Assume with- 
out loss of generality that the vertex at which P starts is and the vertex at which 
it ends is by. There are two things that could go wrong to prevent transforming M 
into MAP from being a shift on M's matching matrix. One is that there exists an 
unmatched spike on either spike train between xy and yy . In that case, it can be 
shown that if M is the minimum matching of cardinality k, then MAP is not the 
minimum matching of cardinality k + 1 ; it can be improved upon by replacing an 
edge. See Figures 5a and 5b. The second possible problem is that there exists an 
edge in M that some edge in P crosses. Then, it can be shown that M or MAP 
contain crossing edges, contradicting Lemma 3. See Figures 6a and 6b. With these 
two possibilities eliminated, P could only be an M-augmenting path corresponding 
to a shift. 

Theorem 3. Let M be a minimum weight matching of cardinality k. Lf a mini- 
mum weight matching of cardinality k+1 has greater weight then M, then M is a 
minimum weight matching over all matchings of any cardinality. 

The proof is in the Appendix. 

Corollary 4. Let denote the minimum weight matching of cardinality k. Lf 
M ko is the minimum weight matching over matchings of all cardinalities, than no 
Mk contains a positive-weighted edge for k < fc . 
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Proof. The proposition is vacuous if kg = 0. So assume the statement is true for k, 
< k < k , we will prove it for k + 1. Assume that M/~ + i has an edge e of positive 
weight. Then Mk+i A{e} is a matching of cardinality k, so w(i4+iA{e}) > w(Mk). 
But w(Affc+i A{e}) < w(Mk+\) < w(Mk) by Theorem 3, a contradiction. □ 

The implication of this corollary is that if any two adjacent columns of a matching 
matrix contain entries that are separated by more than 2 1 / p /q, the matrix can be 
cut in two and the two halves optimized separately. 

7. Conclusion 

We have presented a spike metric satisfying two important desiderata: That it be 
grounded in the time-coding hypothesis of spike generation, and that it be closely 
related to the Euclidean £2 norm. The latter property is important to subsequent 
stages of analysis of spike data, especially if MDS is used. 

We have also presented a fast algorithm to calculate our metric and proved 
its correctness and amortized complexity. We demonstrated that it is faster than 
existing algorithms to compute the p = 1 special case. Our proof used the toolboxes 
of graph theory and combinatorial optimization, and demonstrates that they can 
be usefully brought to bear on problems in computational biology. 

We conjecture that related procedures to the one we have presented could be 
used to align both spike trains with inputs from multiple neurons (as done in 
|Aronov, 2003| ) and strands of DNA. 

Future work will include applications to multiple alignment of neural signals and 
DNA. Our metric suggests new measures of the quality of multiple alignments, and 
our algorithm suggests new ways of computing them. One approach to aligning 
multiple spike trains would be to generalize the CLUSTAL procedures for 'progres- 
sive' alignment, i.e. building a multiple alignment from a collection of aligned pairs 
|Higgins, Sharp, 1988 , |Thompson, Higgins, Gibson, 1994| . 



8. Appendix 

Definitions. G is connected if one can travel from any vertex to any other vertex 
via the edges. A component of G is a maximal connected subgraph. A path in 
which the first and last edges share a vertex is called a circuit. 



Lemma 4. |Schrijver, 2 003 , Theorem 1: If M is a matching in G then either there 



exists an M -augmenting path P or no matching of greater cardinality exists. If N 
is a matching of greater cardinality, we can pick P C M U N . 

Proof. If M is a matching of maximum cardinality and P is an M-augmcnting path 
then MAP would be a matching of greater cardinality, a contradiction. If M is not 
a matching of maximum cardinality, there exists a larger matching N . Let H be 
the graph with vertices V and edges M U N. Every connected component of H is 
a path or a circuit. As \N\ > \M\ one of these components C contains more edges 
of N than M. Its edges must alternate between those in M and those in N, since 
no two edges of either one can be connected to the same vertex. Thus, in order to 
have more edges in TV than in M, it must start and end with non-identical edges 
in N. Therefore, C is an M-augmenting path. 

Proof of Theorem 1. From Schrijvc r, 2003] , Section 3.5, Proposition 1: Let N be 
an arbitrary matching of size k + 1. By Lemma 4 we can pick an M-augmenting 
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path Q in M U N. By definition l(Q) < l(P). \NAQ\ = k, and since M is of 
minimum weight for cardinality k, w(NAQ) > w(M). Therefore, 

w(N) = w(NAQ) - l(Q) > w(M) - l(P) = w(M') 

Therefore M' is of minimum weight over matchings of cardinality k + 1. 

Proof of Lemma 2. Let a = — yj 2 , = Xi 2 — yj 1 , and 7 = yj 2 — yj 1 . Now, 
< 7 < — a. We want 

/(7) = M p + |/3| p - |a + 7l p - 1/3 ~ 7l p > 0. 
The left hand side has roots for 7 at and — a. 

f\l) = -p-sign(a + 7)|a + 7 | p - 1 +p ■ sign(/3 - j)\0 - ^p- 1 . 
This expression is zero if and only if 7 = (0 — ot)/2. By Jensen's Inequality, 

/(^)=M P +l/3| P -2 

The lemma follows by Rolle's Theorem. 

Proof of Theorem 2. We assume without loss of generality that P is of the form 

({ a n . b n }, {a ll ,b j2 },..., {a lt , bj, }). 

Assume for the sake of contradiction that for some t either i t < i t+1 and i t+ i > i t+2 
or i t > it+i and i t +i < it+2 (The cases where for some t either j t < j t +i and 
jt+\ > jt+2 or jt > jt+i and j t +\ < jt+2 arc analagous). The cases are similar so 
we only treat the first one. It has two subcases: jt+i > jt+2 and j t +\ < jt+2- In the 
first subcase, the edges {di t ,bj t+1 } and {<ii t+1 , bj t+2 } are both in M, contradicting 
Lemma 3. In the second subcase case, the edges {dj t+1 , bj t+1 } and {fli t+2 , bj t+2 } are 
both in MAP, which by Theorem 1 is a minimum weight matching of cardinality 
k + 1, contradicting Lemma 3. 

Proof of Main Theorem. We show that if M is a minimum weight matching of 
cardinality k then an M-augmcnting path of maximum length, P, is of the form 

({a», bj}, {a t , bj+i}, . . . , {ai+ N ,b j+N }) or {{a t , bj}, {a i+ i, bj}, . . . , {a i+N , b j+N }) 

where in the first case there are no unmatched elements of x between Xi and yj and 
there are no unmatched elements of y between Xi+M and yi+N, and in the second 
case there are no unmatched elements of y between x% and yj and there are no 
unmatched elements of x between x i+ N and yi+N- Therefore, transforming M into 
MAP is equivalent to the shift operation described in Section 4. 

We first prove the case where p > 1. Using Theorem 2, assume without loss of 
generality that P is of the form 

(Ri , b h L Ri > h 32 }> • • • i i b j, » 
where {ai t } and {bi t } are strictly increasing. The theorem is clear when 1 = 1. So 
let Z > 1 and assume that there exists and integer r such that i t < r < i t +\ (we 
can likewise assume that there exists an integer between j t and jt+i, this case is 
similar). Now there are two cases: either a r is connected to an edge in M or it 
is not. If it is connected to an edge, call that edge {a r ,b s }. s > jt+i, otherwise 
the fact that {a r ,b s } and {ai t ,bj t+1 } are both in M contradicts Lemma 3. But in 
that case, since {a ri b s } ^ P, {a r ,b s } and {di t+1 , bj f+1 } are both in MAP, which 
by Theorem 1 is a minimum weight matching of cardinality k + 1 , contradicting 



a + 



> 
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Lemma 3. 

Now consider the case when a r is not connected to an edge. Replace the edge 
{<2j t+1 , bj t+1 } with {ct r ,bj t+1 } in P to form P' . P' may not be an M-augmenting 
path, but MAP' is still a matching of cardinality k + 1. We know that g p |xi t — 
2/jt+i I P ^ Q^l^r — 2/j't+i I P since otherwise the edge {a r , bj t+1 } would be in M instead 
of {a lt ,b jt+l }. Therefore, \x H - y jt+1 \ < \x r - y H+1 \. If y H+1 > x r , this implies 
—Xi t + Vj t+1 < —x r + Uj t+1 making Xi t > x r , a contradiction. So x r > Uj t+1 - Since 
x it+1 > x r , x r - y Jt+1 < x it+1 - y h+1 , making 

w(MAP') - w(MAP) = qP\x r - y 3t+1 \ p ~ <f K +1 - y 3t+1 \ P < 0- 

MAP' is a matching of cardinality k + 1 of lower weight than that of MAP, but 
since P is of maximum length, Theorem 1 implies that MAP has minimum weight. 
This is a contradiction and the theorem follows for p > 1 . 

Finally, assume that there is an unmatched element of x between x^ and y 3l , 
call it x u , with corresponding vertex a u (the case when there is an unmatched 
element of y between x^ and y 3l is analogous). Since x u < x ix by the previous 
paragraph, yj ± < x u < Xi 1 . Replace the edge {a^^j-^ with {au^j^ in P to form 
P'. P' may be be an M-augmenting path, but as before MAP' is still a matching 
of cardinality k+1. x u < Xi ± implies x u — yj 1 < x^ —yj 1 , and since these quantities 
are nonnegative, 

w(MAP') - w(MAP) = qP\x u - y n \ p - q p \x n - y n f < 0. 

Using Theorem 1, P could not have been an M-augmenting path of maximum 
length, a contradiction. The theorem follows for p > 1. 

Now we address the p = 1 case. For some constants C G Z, eci, . . . , aj > 0, and 
u\,. .. ,uj £ {±1}, the difference in length between two Af-augmenting paths P 
and P' is 

J 

g(p) = l(P') - l(P) = M »< + 2a 

i=l 

g(p) is real analytic and its Taylor series T(jp) converges everywhere. Either g(p) is 
identically or there exist some least integral exponent e such that the coefficient 
of p e in T(jp) is nonzero. If it is positive then there exists some t\ > such that 
f(p) is strictly positive on (1, 1 + ei), likewise if it is negative we can pick ei such 
that g(p) is strictly negative on (1, 1 + ei) by |Nash, 1959| , Theorem 1. We can 
determine such an e\ for every pair of M-augmenting paths P and P' that do not 
have the same length for all p. Let be the minimum of all of these ei's. There 
exists a set of M-augmenting paths Qi, - ■ ■ , Qk such that for p E (1, e 2 ) they all 
have equal length and they are all longer than all other M-augmcnting paths. They 
all must be shifts. By continuity, when p — 1 they still have equal length and are 
at least as long as all other M-augmenting paths. Therefore, MAQi is a matching 
of minimum weight of cardinality k + 1 when p = 1 . 

Definition. Let M be a matching of minimum weight of cardinality k and P be 
an M-augmenting path. We call P an M-shift it it corresponds to a shift in the 
matching matrix of M. 

Proof of Theorem 3. Let M& be the minimum weight matching of cardinality k. 
It suffices to show that if w(Mk) < w(Mk+i), then if M^+2 exists, w(Mk+i) < 
w(Mk + 2)- It follows by induction that w(M} I ) has lower weight than all min- 
imum weight matchings of cardinality greater than k. Therefore, the sequence 



A FAST C p SPIKE ALIGNMENT METRIC 



11 



w(Mi), w{M2), ■ ■ ■ , iw(M m i n ( m „)) is decreasing then increasing, and the point at 
which the switch occurs is its minimum. 

All Mfc-shifts have negative length, since the one with the greatest length, P, 
has negative length, as w(Mk) < w(Mk+i) (using Theorem 1). Therefore all Mk+i- 
shifts have negative length since they are also Mk shifts, except for possibly one, 
Q, that is not. If Q does not exist or is not the Mfc+i-shift of greatest length, 
then the shift of greatest length has negative length, and w(M k+ i) < w(M k+2 ), 
so we are done. So assume that Q is the Mfc+i-shift of greatest length. We shall 
prove that l(Q) < 0. Sending Mk to M^AP = M k +i and then sending Mk+i to 
Mk+iAQ = Mk+2 corresponds to deleting two pairs of stars in one of four classes 
of matching matrices. We consider only two of them, since the others are analogous: 

x i x i+l ' ' ' * ' ' ' ^i+s • • • %i+t * 

_ ■■■ * V] ■■■ V]+r ■■■ * ■■■ Vj+t-1 Vj+t ■■■ _ 

• • • Xi ' ' ' %i-\-r ' ' ' * ' ' ' %i-\-t * ' ' ' 

. ••• * Vj ■■■ * ■■■ Vj+s ■■■ Vj+t-i Vj+t 

To get the other two cases, switch the top stars to the bottom and the bottom 
stars to the top in each of the above examples, and rcindcx accordingly. In both 
matrices, deleting the inner pair of stars corresponds to transforming M k into Mk+i 
and deleting the outer pair of stars corresponds to transforming Mk+i into Mk+2- 
In the first matching matrix there are three possible M^-shifts. The middle one is 
P. Call the rightmost one P_ and the rightmost one P + . l(P+) < l(P) < 0, and 
l(P-) < l(P)- l{P) + l{Q) = l{P-)+l(P+), since performing the shift P followed by 
Q is the same operation as performing the shift P_ followed by P + on the matrix. 
Therefore l(Q) = l(P+) - l(P) + Z(P_) < 0. 

Now consider the second matching matrix. There are M^-augmenting paths P_ 
and P + such that changing Mk to M^AP_ or MfeAP + is equivalent to deleting 
both the first and third stars in the matrix (from the left) or both the second and 
fourth stars from the matrix. As before, Z(P+) < l(P) < 0, l(P-) < 1{P), and 
l(P) + l(Q) = l(P-) + l(P+). Again, l(Q) = l(P+) - l(P) + l(P_) < 0. 

References 

[Aronov, 2003] D. Aronov, "Fast algorithm for the metric-space analysis of simultaneous responses 
of multiple single neurons," Journal of Neuroscience Methods, Vol 124, Issue 2, 
(April, 2003) 175-179. 

[Aronov et al., 2003] D. Aronov, D. Reich, F. Mechler, and J. Victor, "Neural coding of spatial 
phase in VI of the macaque monkey," Journal of Neurophysiology 89: 3304-3327, 
2003. 

[Aronov, Victor, 2004] D. Aronov and J. Victor, "Non-Euclidean properties of spike train metric 
spaces," Physical Review E 69 (2004) 061905. 

[Burkard, 1999] R. Burkard, "Selected topics on assignment problems," Discrete Applied Mathe- 
matics, Volume 123, Issue 1-3 (November 2002): 257-302. 

[Burkard, Klinz, Rudolf, 1995] R. Burkard, B. Klinz, and R. Rudolf, "Perspectives of Monge prop- 
erties in optimization," Discrete Applied Mathematics, Volume 70, Issue 2 (Septem- 
ber 1996): 95-161. 

[Carrillo, Lipman, 1988] H. Carrillo and D. Lipman, "The multiple sequence alignment problem 
in biology," SIAM Journal on Applied Mathematics, Vol. 48, No. 5 (Oct, 1988): 
1073-1082. 



12 ALEXANDER J. DUBBS 12 , BRAD A. SEILER 12 , AND MARCELO O. MAGNASCO 1 



[Chase, Young, 2006] S. Chase and E. Young, "Spike-timing codes enhance the representation of 

multiple simultaneous sound-localization cues in the inferior colliculus," Journal of 

Neuroscience, April 12, 2006: 26(15):3889-3898. 
[Ding, He, 2004] C. Ding and X. He, "K-means Clustering via Principal Component Analysis," 

Proceedings of the International Conference on Machine Learning (ICML 2004), pp 

225-232. July 2004. 

[Gerstner et al., 1997] W. Gerstner, A. Kreiter, H. Markram, and A. Herz, "Neural codes: Firing 
rates and beyond," Proceedings of the National Academy of Sciences, Vol. 94, pp. 
12740-12741, November 1997. 

[Hardle, Simar, 2004] W. Hardle, L. Simar, Applied Multivariate Sta- 

tistical Analysis, Version 26, August 2004, available at: 
www.quantlet.com/mdstat/scripts/mva/htmlbook/mvahtml.html 

[Higgins, Sharp, 1988] D. Higgins and P. Sharp, "CLUSTAL: a package for performing multiple 
sequence algignment on a microcomputer," Gene, 73 (1988): 273-244. 

[Hoffman, 1963] A. Hoffman, "On simple linear programming problems," Proceedings of Symposia 
in Pure Mathematics Volume VII: Convexity, 1963, pp. 317-327. 

[Hopfield, 1995] J. Hopfield, "Pattern recognition computation using action potential timing for 
stimulus representation," Nature, Vol. 376, 6, July 1995, pp. 33-36. 

[Kuhn, 1955] H. Kuhn, "The Hungarian Method for the assignment problem," Naval Research 
Logistics Quarterly, 2:83-97, 1955. 

[Kruskal, Wish, 1978] J. Kruskal, M. Wish, Multidimensional Scaling, Sage University Paper se- 
ries on Quantitative Application in the Social Sciences, 07-011. Beverly Hills and 
London: Sage Publications (1978). 

[Lipman et al., 1989] D. Lipman, S. Altschul, and J. Kececioglu, "A tool for multiple sequence 
alignment," Proceedings of the National Academy of Sciences, Vol. 86, pp. 4412- 
4415, June 1989. 

[Monge, 1781] G. Monge, "Deblai et remblai," Memoires de l'Academie des Sciences, 1781. 
[Munkres, 1957] J. Munkres, "Algorithms for the Assignment and Transportation Problems," 

Journal of the Society of Industrial and Applied Mathematics, 5(l):32-38, 1957 

March. 

[Nash, 1959] S. Nash, "The Higher Derivative Test for Extrema and Points of Inflection," The 
American Mathematical Monthly, Vol. 66, No. 8, (Oct., 1959), pp. 709-713 

[Notredame, 2002] C. Notredame, "Recent progresses in multiple sequence alignment: a survey," 
Pharmacogenomics (2002) 3:131-144. 

[Papadimitriou, Steiglitz, 1998] , Combinatorial Optimization: Algorithms and Complexity, 
Dover Publications, Mineola, New York, 1998. 

[Rieke et al. 1997] F. Rieke, D. Warland, R. de Ruyter van Steveninck, W. Bialek, Spikes: Ex- 
ploring the Neural Code, The MIT Press, Cambridge Massachusetts, 1997. 

[Roweis, Saul, 2000] S. Roweis and L. Saul, "Nonlinear dimensionality reduction by locally linear 
embedding," Science, New Series, Vol. 290, No. 5500, (Dec. 22, 2000): 2323-2326. 

[Samonds, Bonds, 2003] J. Samonds and A. Bonds, "From another angle: differences in cortical 
coding between fine and coarse discrimination of orientation," Journal of Neuro- 
physiology 91: 1193-1202, 2004. 

[Schrauwen, Campenhout, 2007] B. Schrauwen, J. van Campenhout, "Linking non-binned spike 
train kernels to several existing spike train metrics," Neurocomputing, Volume 70, 
Issues 7-9, March 2007, Pages 1247-1253. 

[Schrijvcr, 2003] A. Schrijver, A Course in Combinatorial Optimization, available at 
http: / /www. math. ku.edu/~jmartin/old-courscs/math996-F06/Schrijver.pdf 

[Sellers, 1974] P. Sellers, "On the theory and computation of evolutionary distances," SIAM Jour- 
nal on Applied Mathematics, Vol. 26, No. 4 (Jun, 1974): 787-793. 

[Thompson, Higgins, Gibson, 1994] J. Thompson, D. Higgins, and T. Gibson, "CLUSTAL W: im- 
proving the sensitivity of progressive multiple sequence alignment through sequence 
weighting, position-specific gap penalties and weight matrix choice," Nucleic Acids 
Research, 1994, Vol. 22, No. 22, 4673-4680. 

[van Rossum, 2001] M. van Rossum, "A novel spike distance," Neural Computation 2001 
Apr;13(4):751-63 

[Victor, 2005] J. Victor, "Spike train metrics," Current Opinion in Neurobiology 2005, 15:585-592. 



A FAST C p SPIKE ALIGNMENT METRIC 



13 



[Victor, Purpura, 1996] J. Victor and K. Purpura, "Nature and precision of temporal coding in 
visual cortex: a metric-space analysis," Journal of Neurophysiology, Vol. 76, No. 2, 
August 1996: 1310-1326. 

[Victor, Purpura, 1997] J. Victor and K. Purpura, "Metric-space analysis of spike trains: theory, 
algorithms, and applications," Network 8, 127-164 (1997). 



14 ALEXANDER J. DUBBS 12 , BRAD A. SEILER 12 , AND MARCELO O. MAGNASCO 1 




Figure 1: M-augmenting paths. 

(a) a matching M of cardinality 3 between two sets (red and green) 
representing two spike trains. 

(b) An M-augmenting path P, with edges that already were in M in black and 
new edges in blue; the path should start and end in blue (new) edges, and 
alternate between blue (new) and black (old); it need not contain any black 
edges at all. 

(c) A new matching MAP is obtained as the symmetric difference between the 
first two graphs. 
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Figure 2: Shifts are a subset of M-augmenting paths. 

The matching on top has been drawn with the convention stated below Lemma 1: 
Red nodes corresponding to spikes in train x are on top, green nodes 
corresponding to spikes in train y are on bottom, and nodes farther to the right 
correspond to spikes occurring later in time. All future figures also follow this 
convention. 

The M-augmenting path P on the bottom contains, as in Figure 1, black edges in 
M and new blue edges. Because P moves monotonically forward in time and does 
not skip any node, it is a shift. 
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Figure 3: Convexity forbids edges from crossing. The matching depicted 
on the left cannot be a minimum weight matching of cardinality 2 because 
it violates Lemma 3, since it has crossing edges. The matching depicted on 
the right could be a minimum weight matching of cardinality 2, as it does 
not violate Lemma 3. 




Figure 4. A shift cannot turn back on itself, because either the initial matching or 
the final matching would have crossing edges. 

(a) M is the matching depicted by the two dotted black lines. The complete set of 
lines forms an M-augmenting path P. M is not a minimum weight matching of 
cardinality 2 by Lemma 3 contradicting the hypothesis of Theorem 2, so we do not 
need to consider this case in its proof. 

(b) Here, M is the matching of cardinality 2 depicted by the two black lines, one 
dotted and one solid. The two blue lines and the one dotted black line form an M- 
augmenting path. MAP is the pair of blue lines and the one solid black line. MAP 
contains crossing edges, so by Lemma 3 it is not the minimum weight matching of 
cardinality 3. Thus, by Theorem 1, P is not an M-augmenting path of maximum 
length. 
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Figure 5. Shifts cannot skip previously unmatched vertices. 

(a) : M is the matching composed of the two black lines, and the M-augmenting path 
P comprises those edges plus the three blue ones. P cannot be maximum-length 
because it skips the yellow vertex. The skipped vertex can be joined to a unique 
green vertex without crossing lines; we color it green. The green edge must be 
shorter than one of its two neighbors in the path; if it was shorter than the black 
edge, then the original matching would not be minimal. If it is shorter than the blue 
edge, then replacing the green edge for the blue edge would give a new matching 
that is better than the old one, but the difference is not an M-augmenting path; by 
Theorem 1, there must be an even better one. 

(b) : the last edge in a path cannot jump over a vertex. M is the single dotted black 
edge, and the M-augmenting path P comprises that edge and the two blue ones. 
MAP would have lower weight if the green edge replaced its nearest blue edge in P, 
so P is not the M-augmenting path of maximum length, by Theorem 1. (In this figure, 
the optimal shift would be just the green edge). 
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Figure 6. Shifts cannot skip previously matched vertices, because the edges from 
these previously matched vertices must cross either an edge in the old matching or 
in the new matching. 

(a) M is the matching composed of the three black lines, where the dotted line 
represents an old edge not in the path; and an M-augmenting path is superimposed. 
The dotted line crosses a black line in the path, therefore the old matching was not 
minimal. 

(b) same convention, M is the matching composed of the two solid black lines and 
the dotted black line. An M-augmenting path P is composed of the three blue lines 
and two solid black lines. MAP is the matching composed of the three blue lines and 
the one dotted black line. MAP has crossing edges. 



